template< typename T, unsigned int R, unsigned int C >
Matrix<T,R,C>::Matrix() {
  T* ptr = (T*)elems;
  for( unsigned int i = 0; i < R*C; ++i ) {
    ptr[i] = T();
  }
}

template< typename T, unsigned int R, unsigned int C >
Matrix<T,R,C>::Matrix( T* vals ) {
  T* ptr = (T*)elems;
  for( unsigned int i = 0; i < R*C; ++i ) {
    ptr[i] = vals[i];
  }
}

template< typename T, unsigned int R, unsigned int C >
Matrix<T,R,C>& Matrix<T,R,C>::operator ()( T* vals ) {
  for( unsigned int i = 0; i < N*M; ++i ) {
    elems[i] = vals[i];
  }
  return this;
}

template< typename T, unsigned int R, unsigned int C >
template< unsigned int X, unsigned int Y > 
Matrix<T,R,C>::Matrix( const Matrix<T,X,Y>& matrix ) {
  unsigned int i = 0;
  for( ; i < ((R>X)?X:R); ++i ) {
    elems[i] = matrix[i];
  }
  for( ; i < R; ++i ) {
    elems[i] = Vector<T,C>();
  }
}

template< typename T, unsigned int R, unsigned int C >
Matrix<T,R,C+1> Matrix<T,R,C>::AddColumn( unsigned int target, const Vector<T,R>& val = Vector<T,R>() ) const {
  Matrix<T,R,C+1> temp;
  for( unsigned int i = 0; i < R; ++i ) {
    for( unsigned int k = 0; k < target; ++k ) {
      temp[i][k] = elems[i][k];
    }
    temp[i][target] = val[i];
    for( unsigned int k = target; k < C; ++k ) {
      temp[i][k+1] = elems[i][k];
    }
  }
  return temp;
}

template< typename T, unsigned int R, unsigned int C >
Matrix<T,R,C-1> Matrix<T,R,C>::RemoveColumn( unsigned int target ) const {
  Matrix<T,R,C-1> temp;
  for( unsigned int i = 0; i < R; ++i ) {
    for( unsigned int k = 0; k < target; ++k ) {
      temp[i][k] = elems[i][k];
    }
    for( unsigned int k = target; k < C-1; ++k ) {
      temp[i][k] = elems[i][k+1];
    }
  }
  return temp;
}

template< typename T, unsigned int R, unsigned int C >
Matrix<T,R+1,C> Matrix<T,R,C>::AddRow( unsigned int target, const Vector<T,C>& val = Vector<T,C>() ) const {
  Matrix<T,R+1,C> temp;
  for( unsigned int i = 0; i < target; ++i ) {
    temp[i] = elems[i];
  }
  temp[target] = val;
  for( unsigned int i = target; i < R; ++i ) {
    temp[i+1] = elems[i];
  }
  return temp;
}

template< typename T, unsigned int R, unsigned int C >
Matrix<T,R-1,C> Matrix<T,R,C>::RemoveRow( unsigned int target ) const {
  Matrix<T,R-1,C> temp;
  for( unsigned int i = 0; i < target; ++i ) {
    temp[i] = elems[i];
  }
  for( unsigned int i = target; i < R-1; ++i ) {
    temp[i] = elems[i+1];
  }
  return temp;
}

template< typename T, unsigned int R, unsigned int C >
Matrix<T,C,R> Matrix<T,R,C>::Transpose() const {
  Matrix<T,R,C+1> temp;
  for( unsigned i = 0; i < R; ++i ) {
    for( unsigned k = i+1; k < C; ++k ) {
      temp[i][k] = elems[k][i];
      temp[k][i] = elems[i][k];
    }
  }

  for( unsigned int i = 0; i < ((R>C)?C:R); ++i ) {
    temp[i][i] = elems[i][i];
  }

  return temp;
}

template< typename T, unsigned int R, unsigned int C >
bool Matrix<T,R,C>::operator==( const Matrix<T,R,C>& compare ) const {
  for( unsigned int i = 0; i < R; ++i ) {
    for( unsigned int k = 0; k < C; ++k ) {
      if( elems[i][k] != compare.elems[i][k] )
        return false;
    }
  }
  return true;
}

template< typename T, unsigned int R, unsigned int C >
bool Matrix<T,R,C>::operator!=( const Matrix<T,R,C>& compare ) const {
  return !( *this == compare );
}

template< typename T, unsigned int R, unsigned int C >
bool Matrix<T,R,C>::IsNull() const {
  for( unsigned int i = 0; i < R; ++i ) {
    for( unsigned int k = 0; k < C; ++k ) {
      if( elems[i][k] != T() )
        return false;
    }
  }
  return true;
}

template< typename T, unsigned int R, unsigned int C >
Vector<T,C>& Matrix<T,R,C>::operator[]( unsigned int target ) {
  return elems[target];
}

template< typename T, unsigned int R, unsigned int C >
Vector<T,C>  Matrix<T,R,C>::operator[]( unsigned int target ) const {
  return elems[target];
}
template< typename T, unsigned int R, unsigned int C >
Matrix<T,R,C> operator-( Matrix<T,R,C>& a) {
  Matrix<T,R,C> temp(a);
  for( unsigned int i = 0; i < R; ++i ) {
    temp[i] = -temp[i];
  }
  return temp;
}

template< typename T, unsigned int N >
T Determinant( Matrix<T,N,N>& a) {
  T determinant = 0;
  Matrix<T,N-1,N> shrunk = a.RemoveRow(0);
  for( unsigned int j = 0; j < N; ++j ) {
    Matrix<T,N-1,N-1> temp = shrunk.RemoveColumn(j);
    determinant += ((j%2)? -1 : 1) * Determinant(temp) * a[j][0];
  }
  return determinant;
}

template< typename T >
T Determinant( Matrix<T,2,2>& a) {
  T* source = (T*)(&a[0]);
  return source[0]*source[3] - source[1]*source[2];
}